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(Abstract) This paper revisits the problem of a spinning top in a uniform gravitational field when one point on the symmetry axis is 
fixed in space. It is an instructive and synthetic work of which the theoretical part includes all necessary issues to formulate the full 
differential equations governing the general motion of the spinning top under arbitrary initial conditions. Both Euler and Lagrange 
formulations are discussed. Moreover, closed form analytical solutions are derived for the regular precession and the nutation. The 
numerical integration of the equations was achieved using the standard Runge-Kutta scheme ODE45 available in MATLAB®, which 
was initially applied to the totality of Euler' s equations and then to Lagrange's equations. Also, in house RK2 and RK4 Runge-Kutta 
as well as Crank-Nicolson schemes were applied in conjunction with the constraint for energy conservation. The quality of the 
numerical solution was evaluated by testing the conservation of total energy as well as angular momenta in the form of residuals in the 
corresponding Euler' s dynamic equations. 
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1. INTRODUCTION 

Symmetric spinning top is a case studied by Lagrange in 1788 
[1]. After 1 10 years, in 1897-1898 Klein and Sommerfeld [2,3] 
published two volumes, whereas in 1897 Klein separately 
published a shorter monograph on the mathematical theory of 
the top [4]. In the eastern world literature, one of the oldest 
papers is probably due to Appel'rot ( 1 894) [5] . Historical books 
of general interest are [6,7], whereas another book by Gould [8] 
includes 367 references until early 1970s. Explicit integration 
of motion equations, to give the nutation in terms of elliptic 
integrals is cited in the aforementioned books at the end of 
nineteenth century [3,4] as well as in Whittaker [9]; some 
useful explanations are due to Zaroodny [10]. 

Older standard books of physics and classical dynamics 
[1 1-14] include analytical formulas for nutation, whereas recent 
physics standard textbooks [15,16] limit their discussion only 
to the case of regular precession. There is no doubt that the 
spinning top, as a toy, keeps its fascination from little children 
to adults. A simple internet search reveals many references 
concerning the spinning tops and gyros as if they were 
"magical" instruments that defy gravity. Even the renowned 
mathematician Michele Audin has used the term "mysterious 
(?)" when referring to Kowalevski case in the introduction of 
her book [17]. It is worth-mentioning that, in 1889, 
Kovalevskaya showed that the rigid body motion was 
integrable under certain conditions concerning the ratio (1 :2) of 
the principal momenta between other parameters; her work was 
so remarkable that it won her the Bordin prize (1888) 
[12,34-37]. Besides the toy, there are many industrial 



applications such as navigation of the closely related gyroscope. 
However, a detailed literature survey on the gyroscope and its 
applications is outside the scope of this paper. 

The integrability of Euler's equations describing the motion 
of a spinning top has become a matter of intensive research 
within the last fifteen years. In brief, the numerical solution 
may sometimes not fulfill the law of energy conservation or 
may suffer from "gimbal lock" singularities concerning the 
Euler angles [17-19]. In the regime of numerical analysis, the 
first paper in the western literature referring (among others) to 
the numerical integration of the differential equations of a 
spinning top is probably due to Gorn [20,p.79]. Later, McGill 
and Long [21] studied the case of an unsymmetrical rigid body. 
Simo and associates have developed numerical schemes to 
preserve energy and momentum [22,23] and papers therein. 
Ratiu and Moerbeke [24] have discussed the same matter with 
main emphasis put on the symplectic structure of the motion. 
Historical references have been recently given by Romano [25], 
whereas most recent publications are [26-31]. 

Despite the abovementioned progress, the applicability of 
general purpose numerical integration schemes remains an 
open issue. This paper contributes in this direction by 
developing two variations of the differential equations and then 
implementing standard MATLAB functions such as ODE45, 
which is a code based on a pair of one-step explicit 
Runge-Kutta formulas. The study investigates the performance 
of this time-integrator for the conservation of the energy and 
the angular momentum. For the sake of briefness, the study 
reduces to the symmetrical spinning top. 

2. EQUATIONS OF MOTION 



2.1. Euler's equations 



The top is a rigid body fixed at the point O. Its position 
(orientation) at any instant can be fully described by three Euler 
angles. Following Targ [13], the fixed co-ordinate system 
(space axes) is denoted by (xi^y.^.i) whereas the rotating 
system (body axes) is denoted by (x,y^). The unit vectors for 

the space and body systems are denoted by (ipj^kj and 

(i,j,k) , respectively. According to the original definition, 

first we define the line of nodes (ON) as the intersection 
between the Oxijj and the Oxy coordinate planes; in other 
words, the line of nodes is the line perpendicular to both z.\. and 
z axis and therefore it always lies on the fixed horizontal plane 
Oxijj. The unit vector along (ON) is denoted by n and it is 
chosen so as the system (k) ,k,n ) is right-handed. As usual, we 
define the Euler angles as: 

• is the azimuthal angle between the xi -axis and the line 
of nodes (ON). 

• is the inclination (lean) angle between the zi-axis and 
the z-axis. 

• y/ is the spin angle between the line of nodes (ON) and 
the x-axis. 

Different authors may use different sets of angles to describe 
these orientations, or different names for the same angles, 
leading to different conventions (see, for example, Berry and 
Shukla [32] or the textbook by Hay [33]). 

Let C be the center of mass of the spinning top, which is 
taken along the z-axis at a distance / from the fixed origin O. 
The two first Euler angles ( <f> and ) are directly related to the 

polar angles <f) p and p (known as two of the usual spherical 

coordinates: (f) p ,0 p ,r c — I) of the axis, whose direction is the 

unit vector ki, as follows: (j) =(/) — 7t/2 and p — 0. The 
angle y/ describes the rotation of a material rotation of a 
material point P about the axis of the top, measured relative to 
the intersection of the top with the instantaneous constant 

<f) p -plane. The latter is better understood if at time t = we 

assume the body x-axis to coincide with the line of nodes (ON); 
at this case the plane (Oyz) is perpendicular to the plane (Ox\y\) 
and it intersects the circular cross section at the point P 
(Figure 1). Irrelevant to a possible nutation, it is obvious that at 
a later instant t, the material point P will have rotated along the 
circle (transverse to the body z axis) exactly by the Euler angle 

¥ ■ 

Since the spinning top has a fixed point at O, at every instant 
the motion of the rigid body is a pure rotation about the 
instantaneous vector of the angular velocity CO (along a line 
OQ, not shown), which can be expressed in terms of the Euler 

angles as follows: t» = (/> , C0 e —0, CO^ = ys ■ In more details, 

CO. is along the space z\ axis, CO g is along the line of nodes 

(ON), while 00^ is along the body z axis. As known [1; 12 
p. 174], adding these components of the separate angular 
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velocities, the components [co x ,co y ,(oA of CO with respect to 
the rotating body axes Oxyz are given by: 

co r =0sin#sini// + 9 cosy/ , 

co =0sin#cosi//-#sini//, (1) 

co_ = (f> cos 9 + y) . 

If now the body axes (xyz) are taken as the principle axes 
(123) relative to the reference point O, with moments of inertia 
(Zi, 1.2, h), the abovementioned components of ra are now 

denoted by (CO^CO^CO^ ). Then, the totality of Euler equations 

(kinematic and dynamical) is given as follows: 

i) Kinematic Euler equations 

co, =<f> sin sin y/ + cos y/ , (2a) 
co 2 -<f> sin cos y/ - sin yr , (2b) 
co } -(f) cos + y/. (2c) 

ii) Dynamical Euler equations 

Zjfflj -(/ 2 -7 3 )o 2 co 3 =M l , (3a) 
I 2 cb 2 -(/ 3 -I 1 )co 3 co 1 =M 1 , (3b) 
Zjdjj - (ij - 1 2 ) co 1 co 2 = M } . (3c) 

A rigorous proof of Equations (3) is given in Appendix A. 
The transformation of the angular velocities from the body to 
the space (fixed) system is given in Appendix B. 



*z 
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Figure 1 . Euler angles of a spinning top. 



Equations (2) and (3) provide the general mathematical 
formulation for the motion of the rigid body having one fixed 
point. In the general case of unsymmetrical shape, their solution 
constitutes a difficult mathematical problem [13]. 

In the particular case of a symmetrical spinning top, it is 
convenient to choose the z-axis to be the axis of symmetry, so 

as /, = /, ^ I } . In this case the solution becomes easier. In fact, 

the position vector for the center of gravity C is given as 
follows: 

Space system (Oxi^iZi): 

r ic = l[sin<J)sm6 -cos0sin# cos#] r (4a) 
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r c = [0 if (4b) 
Therefore, the components of the moment (M u M 2 , M 3 ) in 
the right-hand-side of Eq(3abc) are produced by the moments 
of the weight with respect to the origin O [cross product 
r c x (-mg ) ], which is then projected onto the body x-, y-, and 

z-axes, respectively. In the particular case of a spinning top in a 
constant gravitational field, with respect to the world system 
the weight mg (acting at the center C) is the vector 

mg[0 — l] . Considering the transformation matrix A 

obtained after the abovementioned three rotations about the z i , 
ON andz axes, respectively [12, p. 153]: 



Body system (Oxyz): 



' cos (j> cos y/ — sin (j> cos sin y/ 
-cos^sini// - sin cos 9 cos y/ 

V 



sin <j> sin 



sin^cosi// + cos cos ^ sin i// sin#sini//^ 
-sin^sini// +cos0cos0cosi// sin#cosi// 
-cos0sin# cos# 



the projections of the weight onto the axes of the body system 

are given by the matrix product mg-A-[0 — l] , 
whence the torques appearing in (3) are expressed by 



Remarks 



(6) 



'm; 




sin^cosi//" 




• =mgl- 


-sin#sini// > 








^ J 



1) We notice that the Euler angle y/ appears in both the 
kinetic and dynamical expressions (2) and (3) as well as in 
(6). 

2) Since M, =0, by virtue of the condition /, =I 2 Eq(3c) 

implies that C0 3 = , which means that the component <y 3 
is a constant. 

3) A torque with components M } or M 2 will cause both (D, 

and Q) 2 to change without affecting (O i . 

4) Whatever the external loading is, the initial conditions at 
the instant t = are the following six quantities: 

(0o^o)'( 6 'o^o) and (n^o)' ( 7 ) 
from which the angles ^ and y/ are needed only for 

reference (they do not affect the visual orientation of the 
symmetrical spinning top). In other words, the mechanical 
behavior is determined by the four initial conditions: 



(5) 



(J) a ,9 Q ,9 , i/^o )' wn i cn determine the energy (see below 
Eq(18)). 

5) The measure of the vector of the angular velocity 

CO ( co l , a>, , co } ) is given by: 

II II / 9 2 2\'/ 2 

||co|| = I co[ + a> 2 + ft> 3 J 

= (</> 2 + 6 2 +y/ 2 + 2</)6cos6f 2 

6) The angle a between the abovementioned vector to 

and the body z axis is determined by the direction 

cosine ft> 3 / Co | , which by virtue of Eqs(2c) and (8) 
becomes: 



(8) 



COS cc 



UcosG+y/) 



(<j> 2 +e 2 +y/ 2 + 2(j>e cos e ) 



1/2 



(9) 



2.2. Lagrange's equations 

In terms of the instantaneous angular velocity ro of the rigid 
body that moves about the fixed point O, the kinetic energy can 

be computed according to the formula T — l/H^co 2 , where 

Iq^ is the variable moment of inertia with respect to the 
instantaneous axis of rotation OQ. In terms of the components 
of the angular velocity co (ft), , a> 2 , ft> 3 ) in the body system, the 
kinetic energy is given according to the formula 

T = - [i^co; + I 2 co 2 2 + 7 3 ft) 3 ) (10) 



Substituting (2) into (10) and considering a symmetrical top, 
i.e. 7, —I 2 , the latter takes the form: 



T = ^( 2 +<j> 2 sin 2 0) + ^ at 

Also, the potential energy is given by 

V = mgl cos 6 



(11) 



(12) 



Following Targ [13, pp. 518-519], the Lagrange's equations 
are written as follows: 



a 



d 


f dT^ 


dT 


dt 


K d<f>) 


d(j) ' 


d 


f dT\ 


dT 


dt 


Kde) 


86 ' 


d 




dT 


dt 




dy/ 



(13) 



Q 



and take the equivalent form: 



8L 

— = M 
dq 



Alternatively, we can consider the Lagrangian L — T — V 

d_ 
dt 

Substituting (11) into (13) and assuming that the only load is 
due to the weight W=mg of the spinning top, the right-hand 

sides ( , Q e , Q ) of Eq(13) correspond to the moments 

M z =0, M 0N = Wlsm9 (ON= line of nodes) and M z = 0, 
respectively. As a result, Eq(13) becomes: 

— (/^sin 2 6 +I 3 co 3 cos&) = , (14a) 
dt 

Ifi-lJ) 1 sin0cos0 + / 3 ft> 3 <^sin# = mglsinO (14b) 



7 3 ^ = 
dt 



(14c) 



One can notice that Eq(14a) and Eq(14b) do not include the 
Euler angle if/ . 

It is worth-mentioning that (14a) dictates that the value of 
the component of the angular momentum towards the space 
Z] -axis is a constant {first invariant of the system), i.e 



p = dLjdfy = 1$ sin 2 + 7 3 o 3 cos# = c, , 



(15a) 



where the constant value c.\ in (15a) is directly determined 
in terms of three initial conditions (<f) ,0 ,y/ ), mentioned in 
(7), as follows: 
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c, = 7j0 o sin 2 6 + 1 } (J) Q cos 6 + y/ ^ cos 6 , (15b) 

Equation (14c) depicts that tt> 3 is a constant, which can be 

calculated in terms of three initial conditions (<f> ,0 ,y/ ) using 
the kinematic Euler equation (2c), i.e. 



Q) } =<f> cos# +y/ 



(16) 



Obviously, for any instant t > the component (0 3 is given 

by Eq(2c). Equation (2c) depicts that the aforementioned 
material point P undertakes a relative rotation around the 

rotating z-axis and it also undertakes the transport rotation ^k, 
(namely its projection (j>COs6 along the z-axis). 

Since C9 3 is a constant, the second invariant of the system 
(angular momentum towards the bodyz axis) is: 

8L 



P = = Leo, 



(17) 



Unlike Euler' s equations, these equations define the 

motion only of a symmetrical body for which I X =I 2 . 

However, they are simpler than the totality of Euler' s dynamic 
and kinematic equations [ 1 3 , p . 5 1 9] . 

Usually, a criterion to test the accuracy of the numerical 
integration is the energy conservation, which is given by the 
form (E = T+ V)as follows: 



1 / *2 "2 2 \ 1 2 

E = —I\6 +<p sin 6 \ + —I } a> } + mgl cosO = c 2 



(18) 



As previously, the value c 2 in (18) is directly determined in 
terms of the initial conditions (7) and particularly it requires 

four of them (<j> ,8 ,8 ,y/ ). 

3. A CRITICAL REVIEW ON THE MOTION 
OF THE SPINNING TOP 

In the general motion of the spinning top the inclination (lean) 
angle varies in time. The latter motion is defined as 
"nutation". However, the case of "no nutation", which is also 
called "regular precession", appears a particular interest. At this 
point it is worth-mentioning that given the six initial conditions 
[Eq(7)] the time history in the orientation of the spinning top 
can be calculated in an deterministic way. It was previously 
mentioned that two of the initial conditions, that is the initial 

angles ^ and l/^o do not play any significant role since they 

are only reference values. In other words, the history in the 
position of the spinning top depends on the four initial 

conditions <f> ,0 ,0 ,y/ . The aforementioned four initial 

conditions determine the three system invariants, which are the 
component of the angular momentum along the z.\ axis 



( p = C, ), the component of the angular momentum along the 

z axis ( p = I 3 a> } ) and the total energy (E = c 2 ). 
3.1. Constant inclination angle 

When the inclination angle is constant, we then refer to the 
abovementioned "regular" or "steady" or "smooth" precession. 

The investigation of motions under gravity in which the axis 
of the top makes a constant angle with the vertical has been 
reported for example by Hay [33, pp. 91-94]. However, since 
Hay [33] uses a different body-fixed coordinate system, it 
makes sense to adapt his approach to our traditional xyz body 
system. 

The assumption of a constant angle , whence = , 
causes Eq(2a), (2b) and (2c) to give: 

(D x =§ sin sin y/ , co 2 =<f> sin cos y/ , a> } =<f> cos + y/ . 

(19) 

For this particular case, we designate p — <f> for the 
precession, and S =1// for the spin. 

Multiplying (3a) by COS y/ and (3b) by shll// , then 
subtracting by parts and assuming that (I t = l 2 ), one obtains: 

[(/, ) cosOp 2 + I 3 sp - mgl] sin = 0, (20a) 

Substituting (20a) into the abovementioned (3b) already 
having been multiplied by SUll//, one obtains: 

(/,/>cos^)sin0 =0 (20b) 
Lastly, Eq(3c) becomes: 

/>cos0 + i = O. (20c) 

Equations (20a), (20b) and (20c) are identical with those 
obtained by Hay [33] in a different way. 

Again, the totality of Equation (3) is given by equations 
(20). Below, we distinguish two cases. 

1) One solution of Eq(20) is — . In this case the axis of the 

top is vertical, and the top is said to be "sleeping". Then, 
(20c) becomes p + s = , which means that the precession 
is equal and opposite to the spin. 

2) If is not equal to zero, then sin0 can be eliminated in 

(20a,b). Therefore, (20b) yields p = whereas (20c) then 
yields s — .In other words, p = constant, s = constant. 
Equation (20a) is a relation among the three constants p, s 
and . Hence it appears that we may assign arbitrarily 
values for two of these constants and there will exist a 
corresponding motion of the top with a constant. 
Following Hay [33,p.94], Eq(20a) can be solved in s: 

mgl (/,-/,)/? cos 



s = - 



hP 



(21a) 



Equation (21a) offers the dependency of the spin (s) of the 
precession (p) and the inclination angle ( ). In particular, the 
precession appears in both the denominator and the nominator. 
If the precession is small (p«), we note from (18a) that the 
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second term vanishes whereas the first one becomes large and it 
is approximated by the value 

mgl 



s = ■ 



(21b) 



Moreover, whatever the value of the precession is, if the long 
z axis of the top is horizontal (0=^/2, i.e. COS0 =0 ), 
equation (21b) is not an approximation but an exact solution 
and, if solved in p we obtain: 

mgl 



For = 7r/2: p 



(21c) 



which is the well-known from textbooks of physics [15,16]. 
While Eq(2 1 c) is valid for a horizontal top under any conditions, 
for any other inclination it may hold under the abovementioned 
conditions. 

Although Eq(20a) is a quadratic equation in p, it is not 
offered to determine the precession angular velocity p as is, 
because s depends on p according to Eq(2c). Therefore, 
substituting Eq(2c) into Eq(20a) one finally obtains: 

(l l cos9)p 2 -(l 3 co 3 )p + mgl = 0, e±n/2 (22) 

Equation (22) constitutes the sufficient and necessary 
condition to achieve regular precession (no nutation). It is 
mentioned that Eq(22) can be also obtained from Eq(14b) 

putting 9=0 and eliminating sin in both left and right 
parts. 

In more details, when ^ njl , the quadratic equation (22) 
is valid and its discreminant is 



D = (/ 3 fi> 3 ) 2 - 4IJmg cos 



(23) 



As a real solution requires D > , a "conditio sine qua non" 
for the axial angular velocity tt> 3 to achieve an inclination 
angle is: 

2 4I.lmg cos 9 
<y;>^-^ (24) 

Provided Eq(24) is valid ( D > ), the two roots of (22) are 
given by 



1,(0, 



I 3 co 3 



21, cos 



\ 2/j cos j 



Img 
I, cos 6? 



(25a) 



Therefore, the plus and minus sign in (25a) produce two 
solutions: the so-called "fast" precession and the "slow" 
precession solutions. 

In more details, equation (25a) can be rewritten as follows: 
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27, cos6» v ; 



where 



x = ■ 



4lmgl l cos # 



(25b) 



(25c) 



(7 3 fi) 3 ) 2 

Under certain circumstances, Eq(25b) can be further 
simplified. In fact, if it is not only x < 1 but the rotational 
speed of the spinning top is very large ( x <SC 1 ), i.e. 



I 3 G) 3 » ■ s f4lmgl 1 cos# , (25d) 
the root in (25b) is approximated by 

y/l-x «1 x (25e) 

2 

In such a case, (25b) implies that the fast and slow precession 
is approximated by: 



Fast precession: p 



fast 



Slow precession: p , ~ 



/, cos 6 
mgl 



(25f) 
(25g) 



Remarks 

1) One can notice that while the slow precession does not 
depend on the inclination of the axis of the top (lean angle 
) the same does not hold for the fast precession. 

2) The slow precession is the one most commonly used 
experimentally. "In principle, however, one can also obtain 
a fast precession, although it is difficult to start the 
gyroscope off with just the right motion to achieve it. In 
practice, the possibility of the fast precession can be 
ignored" [11, p. 683]. 

3) In the case of slow precession (25g), the value of p 
corresponding to the negative sign of the radicand in Eq(25a) 
is very small. Therefore the measure of the angular velocity 
co = ||io|| is very close to either of m } and y/ (cf. Eq(2c) 

and Eq(9): COS CC Therefore, the direction of the 

vector CO is very close to the axis of symmetry k of the 
ellipsoid (bodyz axis). 

3.2. Nutzation 

Nutation is the case according to which the inclination angle 

varies in time (6^0). 

At this point, we closely follow Klein and Sommerfeld [3, 
p.222] and Goldstein et al. [12,p.213]. 

Using the two invariant angular momenta ( p^ and p 

towards space z ] and body z axis, respectively), the total energy 
can be written as follows: 



p 2 ( Pj,- P cos O) l 

E = ^+ K * w '-+-L9 2 + mglcos9=c 2 (26) 

2L 2L sin' 9 2 



Solving (26) in , then using the transformation 
u — cos , it can be easily verified that the critical (extreme) 

values of the angle at which the derivative vanishes ( 6 = : 
"turning angles") is given by the roots of the following cubic 
polynomial: 



with 



/ '(w) = au 3 + bu 2 + cu + d = , 
a = 2I x mgl 

b = -co](ll-I l I i )-2I i c 2 
c = 2(c/ 3 ffl 3 -I t mgl) 



(27a) 



(27b) 



d = 2I, 



V 2 



-c. 



Dividing by a , (27a) can be also written as follows: 



where 



u +a l u 2 + a 2 u + a 3 -0, 



-b/a, a 2 =c/a, a i =dja. 



(27c) 



(27d) 



Equation (27c) is easily solved by introducing the auxiliary 
variables: 



3a, - a, 



R = 



9a,GL -27a, -2a, 



54 



(27e) 



The solution is determined by the sign of the discreminant 

D' = Q 3 +R 2 . (27f) 
If D' < , the roots are given by: 



u l = 2yJ-Q cos -p 

V 3 J 



u = 2J-Q cos 



u = 2J-Q cos 



-P + -71 

V 3 3 , 

-J0 + -7T 

U 3 ; 



--a. 



— a, 

3 



(27g) 



with 



(27h) 



Therefore, we obtain one or two roots between the interval 
[-1,1] and another non-physical root (|«|>l) . The two 



physically meaningful roots U V U 2 will yield the span of angles 

of the nutation ( < COS 6 <u 2 ). Obviously, when the roots 

U { ,U 2 are coincident, there is no nutation. 

The above analysis can be further simplified if at time t = 
we consider the following initial conditions: 9 =9 a , and 

O — <j> — . In this case, the first solution of Eq(27a) is the 

initial point (m, = cos ) , whereas it can be easilyproven that 
for the next solution m 2 Eq(27a) becomes: 



f(u) = (u 1 -u) 



2 ( 2\ (^ + ^ +C l 2 )/ 

al x [\-u Jh \u x -u) 



(271) 

Setting x = «, - u and x x =u l —u 2 , the quadratic equation 

within the brackets in the right hand side of Eq(27i) takes the 
form: 



xj 2 + px x -q = , 



where 



p = p 2 1 a - 2 cos 6> and g = sin 2 9 , 



(27k) 



(270 



Thus the realistic solution is one root of the quadratic 
polynomial as follows: 



(27m) 



Following Goldstein [12, pp. 215-217], if the kinetic energy 
of rotation about the body z axis is very large compared to the 
maximum change in potential energy: 



1 



-I 3 G>3 » 2mgl , 



(28) 



we speak of the top as being a "fast top". With this assumption 
we can obtain expressions for the extent of the nutation, the 
nutation frequency, and the average of the frequency of 
precession. 



Since 



gives: 



2 



2mgl 



, the condition (28) for the fast top 



^»2 
a 



(29a) 



Therefore, except in the case that / «; / , it holds p » q , 
and therefore the solution of eq(27m) is approximated by: 



9_ 
P 
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Neglecting the term 2cos0 o compared to jc 2 la , Eq(29b) 
can approximate the extent of nutation by the simple formula: 

/, 2mgl 



2 sin 6 , 



(29c) 



where is the initial condition. Thus, the higher the angular 

velocity tt> 3 is the smaller the nutation. 

Since the amount of nutation is small, the term — u 2 ) in 

Eq(27i) can be replaced by its initial value, sin 2 6 . By virtue 
of Eq(29c), and considering the variable transformation 
y = x- xjl for / (u ) = , Eq(27i) finally becomes: 



-co, 



y 



(29d) 



From elementary mechanics Eq(29d) dictates that, for the 
fast top, the frequency of the nutation, in terms of its angular 
velocity k, is given by: 



k = - 



-co, 



(30) 



and therefore it increases the faster the top is spun initially. 
Based on the calculated extent Au = u 2 -u l = cos 6 2 - COS 6 l , 
we can determine the extent of the inclination angle 
A0 = 9 2 —0 1 . Then, the variation of the inclination angle is 
approximated by: 

A0 



e(t) = e l + — (l-cosfc). 



(31) 



According to French [11, p. 694], "if the initial conditions are 
varied, different types of nutational motion may occur, but they 
are all understandable in terms of the principles underlying the 
above analysis". 

Due to the nutation, the angular velocity of precession is 
harmonic of the form: 

mgl 



(j) = (l-cos&?) 



I 3 w 3 



(32a) 



of which the amplitude is recognized as the p slow [Eq(25g)]. 

Although the rate of precession varies harmonically with time, 
with the same frequency as the nutation, the average precession 
frequency is: 

<t>=-^ = Ps^ (32b) 
I 3 0) 3 



4. THREE MAIN ALTERNATIVE 
FORMULATIONS AND THEIR 
NUMERICAL INTEGRATION 

4.1. Based on the totality of dynamic Euler's 
equations 

Since the right-hand side of (3) includes terms in Euler angles, 
it becomes necessary to solve in these angles (primary variables) 
instead of the angular velocities. It will be shown that the 
system of equations takes the standard form 



b x = mgl sin y 3 cos y 5 - a u y 2 y b - a ]4 y 2 y 4 

b 2 = -mgl sin y 3 sin y 5 - a 23 y 2 y b - a 24 y 2 y 4 
-a 25 y 4 y b -a 26 y 2 -a 2J y 4 
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(38) 



with 



y = f(t,y), 



(33) 



where y denotes the vector of five variables. 
Actually, substituting (2) and (6) into (3) we obtain a 
second-order system of three differential equations in the three 

Euler angles ( (j), 9, y/ ). Although the standard Runge-Kutta 
procedure would normally lead to six equations of first-order, 
however, since M i = and due to the symmetry =I 2 ), 

Eq.(3c) implies ri) 3 = . Therefore, as Eq(2c) suggests, a 
linear dependency between the rotation about the z-axis and the 
rest two Euler angles: y/ — a> } - <f> cos 9 . 

Therefore, by choosing the auxiliary variables as: 



y 6 = ( °j,-y2 cos y 3 

a l3 = /j cos y 5 sin y 3 
a l4 =/, smy 5 cosy 3 
Oj 5 = -/, sin y. 
a 



(39) 



K 5 



a, = a> 3 (l 2 -I 3 )siny 5 
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and also: 



y l = (/>,y 2 =<j>,y 3 =6,y 4 =0,y 5 =y/ , 
the final system to be solved becomes: 



(34) 



yi 




y 2 


y 2 




r 2 


- h 


' = < 


y* 


y A 




r 4 






co 3 -y 2 cosy, 



(35) 



a 2J = -I 2 sin y 5 sin y 3 
a 24 = I 2 cos y s cos y 3 

a 25 =-I 2 cosy 5 (40) 
a 26 = -a) } (l 3 - /, ) sin y 5 sin y 3 
a 21 = —o) 3 [I 3 —I l ) cos y 5 

It is noted that in the system (34) the Euler angle y/ is an 
integral part because it is involved in (37) and (38). 

4.2. Based on Lagrange's equations 

It will be shown that the system of equations takes the standard 
form y = f(t,y) , where now y denotes the vector of four 
variables. 

If Eq(15a) is solved in (/> it gives: 

(c, -I 3 G> 3 COS0) 



/, sin 9 



(41) 



The functions r.% and r<± that appear in the right-hand side of 
(35) are derived from the solution of the linear system: 



Substituting (41) into Eq(14b) one receives: 
(c, -I 3 co 3 cos 6?) 



a,, 





r, 








r 4 




b, 



sin cos 



(36) 



/, sin 2 6 

(cj -I 3 a> 3 cos 9) 



(42) 



mgl 

sin 6 H sin 6 



In (36) the four matrix elements are given by: 
a u = I j sin j/ 5 siny 3 

a n =I 1 cosy s 

a 2l = I 2 cos y 5 sin y 3 

a 22 = -I 2 sin y s 

Also, the right-hand sides of (36) are given by: 



(37) 



/, sin 2 

As previously, the last equation of the system is: 

if/ = 6) 3 - <f> cos 9 , 
which when combined with (41) gives the following: 
(c, - I 3 co 3 COS0) 



If/ = (l) 3 



/, sin 9 



cos 6 , 



(43) 



(44) 



Therefore, the vector of the four variable is chosen as 

y = bpj 2 '3 ; 3'3 ; 4] r ' where: 



y x =<f>, y 2 = 9, y 3 =6, y 4 =y/ 



(45) 



The system j = f(t,y) is produced when considering (i) 
Eq(41) for y, , (ii) y 2 =y 3 , (iii) Eq(42) for y 3 , and (iv) Eq(44) 
for y A . 
Remarks 

1) Since the Lagrange equations (14) do not directly include 
the Euler angle yj (it appears indirectly in 14c), the 
position of the center of mass and the estimation of the total 
energy can be found in terms of four variables only, i.e. 

<f>,<i>,o,d. 

2) Since ^ is calculated from Eq(41), the conservation of the 
angular momentum — c { is taken for granted. 

3) Equation (42) includes only the Euler angle 9 , therefore it 
is an ordinary differential equation of second order that 
could be separately solved. The numerical integration of the 
latter could be performed using any known scheme such as 
Runge-Kutta or Crank-Nicolson algorithms. Afterwards, 

the calculation of (j) through Eq(41) and of yi through 

Eq(43) or Eq(44) will follow. The integration of (j) and y/ 
could be separately performed. 

4.3. Based on two system invariants 

Without details, it has been written in several textbooks that, 
theoretically, one equation in the system that describes the 
motion of the spinning top could be replaced by the energy 
conservation [12,13]. In the same framework, it is well known 
that in one-dimensional problems the conservation of energy 
can be used to derive the derivative of the motion variable 
(displacement or angle) and then analytically integrate the 
equations [15, p. 354 (Wiley- Toppan, 1966)]. 

In our case, substituting Eq(41) in Eq(18), and then solving 

in 9 , one obtains: 




1 



i> 3 



- mgl cos 9 



(c, -/ 3 o 3 cos#) 



/, sin 9 



(46) 

The analytical integration of Eq(46) has been previously 
obtained using elliptic integrals [1-4]. Under certain 
circumstances, Eq(46) can be also integrated in a numerical 
way. One difficulty is that the sign preserves its value only 
between two successive "turning" points of the nutation 
provided the initial derivative is different than zero. In 

particular, when the initial condition is 9 — , the 
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Runge-Kutta method is not applicable. To make this point clear, 
without loss of generality, let us consider the simplest forward 
Euler method (not applied in this paper): 

n+l =6 n +At6(t n ,6 n ) . Since Eq(46) fulfils the initial 
condition 9 — , the aforementioned recurrence formula will 

always vanish ( = 9 2 = . . . = 9 n = 0). The same conclusion 

holds for Crank-Nicolson and Runge-Kutta schemes. But even 
if the initial conditions are different than zero, Eq(46) is 

applicable only until the next "turning point" at which 6 — . 

At the latter point, the numerical solution crashes and always 
leads to the same stable value (exactly as happen when the 

initial point was 9 — ). 

4.3.1 Scenario 1 

Again, the abovementioned scenario of considering a system of 
three equations of first order y = f(t,y ) , i.e. Eq(41) for (j) , 

Eq(46) for 9 , and eq(44) for yr , is applicable for only a small 

period of time (less than one-quarter of the nutation period). 

4.3.2 Scenario 2 

Alternatively, we propose the use of Eq(42), which is a second 
order ODE that can be easily solved using either of the 
Runge-Kutta or Crank-Nicolson methods. For every time step, 

we accept the calculated value 9 n+l "as is" and then we modify 
the value of 9 n+l according to Eq(46). Then, we perform the 

next step based on the corrected value of 9 n+1 , and so on. 

In both above scenarios, due to the calculation of (j) 
according to Eq(41), not only the energy (c 2 ) but also the 
angular momentum (c.i) are preserved. 

5. SUPPORT FORCES 

Applying second Newton's law with respect to the Cartesian 
space system, the three components of the support force at the 
fixed point (axis origin O) are given by: 

(47a) 



F . = mx, „ , 

support,* l,C ' 

F 



support, y 
support, z 

= m(g + z lc ) 



(47b) 
(47c) 



By virtue of Eq(4a), the components of the acceleration at the 
'^center of mass C, in the space axes, are given by: 



x iC = I cos (psinO + O sin <f> cos 6 

-(J) 2 + 9 2 ) sin (j) sin 9 + 2<j>9 cos <j> cos 
y\ c = I [0 sin <f> sin 6 - 6 cos <f> cos 6 

+ {J) 2 + 9 2 ) cos (j) sin 6 + 2<j>6 sin^cos#J 
z\ =-/(0sinfl + fl 2 cos#) (48c) 



(48a) 



,, (48b) 



In addition, the two components on the horizontal plane can 
be split in the radial and circumferential direction (with 

direction cosines e r =[sin^ -COS^] r and 
e e = n = [cos^ sin^] r , respectively) as follows: 



F 



support,/ 



F 



support, (9 



sin^ -cos^ 
cos^ sin^ 



F. 



support,.v 



F 



support, v 



(49) 



Substituting (47) and (48) into (49), one finally obtains: 



F = ml 

A support ,r 1 



cos e-(<i> 2 +9 2 ) sin# 

F su PP ort,e =»»/(# sin + 2<j>0 cos 6 ) , 
g-/(6»sin# + # 2 cos6>) 



F — fti 

support, z 



(50a) 
(50b) 
(50c) 



Although the derivative of the Euler angle y/ and the relevant 

component ftJ> 3 influence the numerical solution, it can be 
noticed that they are not included in the above expressions. 

6. NUMERICAL RESULTS 

The theory will be validated by two numerical examples, the 
first taken from literature whereas the second is of practical 
mechanical engineering importance. In addition to ODE45 
(MATLAB), in house Runge-Kutta and Crank-Nicolson 
methods have been used. 

The accuracy of the numerical solution has been tested either 
by the degree of approximating (i) the conservation of energy 
or (ii) the conservation angular momentum. Since Eq(2c) is the 
standard one that always participates in the system of equations, 
it will be completely fulfilled. Moreover, Equations (2a) and 
(2b) are by-products and there is no apparent measure to 
evaluate their accuracy. Consequently, two apparent residuals 

(7? p i? 2 ) are those appearing in the first two dynamical Euler 

equations [(3a) and (3b)], which are given below: 



where 

ci), = sin 9 sin if/ + 9 cos if/ +<f>6 cos 9 sin if/ 

- 9y/ sin if/ + y/(/> sin 9 cos if/ 

a> 2 =(/> sin 9 cos if/ - 9 sin if/ + <f>9 cos 9 cos if/ 

- 9y/ cos if/ — i//(j) sin 9 sin if/ 

Also, the residuals for the violation in the conservation of space 
zi-momentum (R.3) and the conservation of energy (R4) are 
given by: 



[',< 



'sin 2 9 +/3O3 cos 6?] - 
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(51c) 



R, 



(5 Id) 



— /, 1 9 +<b sin" 9 + —Leo, + mgl cos 9 
2 V ' 2 ' ' 

6.1. Example 1: A slow top 

For reasons of comparison, we closely follow the data 
previously used in literature [30, pp. 138-139]. In addition, we 
also present many comments related to the mechanics involved 
and the auxiliary use of the available closed form analytical 
solutions. 

It is reminded that in matrix form (see Appendix A), the 
relationship between the transformation matrix A and the 
matrix of angular velocities CI is as follows: 








-co 3 


m 2 


A r = A r ■ Q, where £2 = 


CO, 





-co 




-co. 


cw, 






The matrix A includes only the Euler angles whereas the 
angular velocity includes their derivatives. Therefore, the initial 
conditions must include data from both matrices. In more detail, 

the initial data where (J) —0,(j) —0, 6 —7t/l6,6 =0, and 
y/ = 7l/l6,y/ = 1 , whence 





1.0 


0.0 


0.0 










A„ = 


0.0 


cos(^ ) 


sin 0„ ) 


and 








0.0 


-sin ) 


cos 


A 


)_ 













-® 3 ,o 






"0.0 


-1.0 


0.0" 


"0 = 


G) 3fi 





-co 10 




1.0 


0.0 


0.0 




.-^2,0 


co l0 







0.0 


0.0 


0.0 



Following [30], the inertia tensor J, the center of mass r c with 
respect to the body system were given as 



(51a) 




"A 





0" 




"7 





0" 










1 








(51b) 


J = 





h 





~ 8 





7 















h_ 










2 



(kg.m ), r c = 



with / = V3/2m. 

Moreover, with respect to the space (global) system, the 
vector of the external torque is taken to be produced by the 
weight mg of the spinning top at the center of mass C, which for 
the fixed (space) system is: 





mg = in 





-9.81 



(m.s~ ), with m = 1 kg. 



As was previously mentioned, the right-hand side of 
dynamic Euler's equations (3) is given by Eq (6), which 
considers the projection of the vector representing the external 
torque onto the body axes. 

This case appears a much extended nutation. The numerical 
results are valid only if the top is mechanically supported by a 
stand that allows it to dip below the horizontal level. In more 
details, the topology of the stand must be so that it is capable of 
exerting either compressive (above the horizontal) or tensile 
(below the horizontal) support forces. The latter (below the 
horizontal) is crucial because otherwise the top will detach 
from its support point. 

6.1.1 Analytical calculations 

A first observation is that the major initial condition 
6 —7t/\6,6 = is concerned with instantaneous stillness of 
the body z axis, which is equivalent with a turning point of 
nutation (u x =7r/l6). The immediately next turning point is 
accurately calculated by Eq(27g), which provides two feasible 
solutions, i.e. «, ~ 0.9808 (corresponds to 1 1 .25 degrees, being 

the initial angle 9 =nj\6) and U 2 ~ -0.9958 (corresponds to 

turn = 174.8 degrees). Therefore, it is anticipated that the 
spinning top will perform large scale nutation in the interval 
6 e [11.25, 174.8] degrees. 

Concerning the angular velocity of precession, it varies 
between zero (at 6 = 1 1 .25 deg) and a maximum value = 

67.97 rad/s (appearing at the lowest position 6 turn = 174.8 

deg). It is noted that if we try to take the first derivative of (j> in 
Eq(41) and equalize to zero, it leads to a quadratic polynomial 

with negative determinant, a fact that proves that (j> is 
permanently positive until the lowest point thus no other local 
maximum exists in between 6 =11.25 and lum = 174.8 
deg. 

It is also noted that in this example, since 
l/2/ 3 ft>3 =0.125«2mg/ = 19.62, Eq(27m) must be applied 

and taking the positive sign we obtain x x « 1.967 , which 

corresponds to 84.8 degrees below the horizontal, a fact 
verified with the numerical solution presented below. As it will 
be validated in the next Section, the set of simple formulas from 
Eq(29c) to Eq(32b) induce tremendous error and therefore are 
not applicable. 

6.1.2 Numerical computations 

In the sequence we present the numerical solution obtained 
using several schemes: Euler [Eq(35)], Lagrange [Eq(45)] and 
two energy conservation schemes. In order to avoid the use of 
elliptical integrals, as an "exact" solution we have taken the 
ODE45 MATLAB solution using an extremely small tolerance 
of order 10" 14 . For a better evaluation, we present the response 
for three cycles of precession ( A0 = 3><360 degrees of 
azimuthal angle), which here corresponds to almost three 



cycles of nutation. For the data of this example, the invariants 

of the mechanical systems are Cj = 0.245 and C 2 = 8.46. 

In Euler formulation (Section 4.1), the first two residuals (R i 
and i? 2 in Eq(51)) are zero whereas the last two (if 3 and R4) 
vary in time (Figure 2). It can be noticed that the period is 
about 2.22 seconds. MATLAB ODE45 was applied using the 
default accuracy. 
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Figure 2. Calculated residuals (R [,R 2 ,Ri,R4) based on Eq(51), using 
ODE45 for the Euler formulation (336 unequal steps). 

In Lagrange formulation (Section 4.2), as previously, the 
first two residuals are zero. In addition, no variation of c \ was 
observed (i? 3 =0). Moreover, not only the energy is not 
preserved (R4 ^ 0) but it also appears singularities when the 
center of mass passes through the horizontal. Concerning the 
angular momentum, both R\ and vary in time (Figure 3). 
MATLAB ODE45 was applied using the default accuracy. 
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Figure 3. Calculated residuals (R i,R 2, i?3,i?4) based on Eq(51), using 
ODE45 for the Lagrange formulation (248 unequal steps). 

In house standard RK2 and RK4 Runge-Kutta algorithms 
using a constant step (Appendix C) during the whole procedure 
leads to extremely high singularities (R] and R2 of the order 
10" 10 , whereas R4 of the order 10 9 ) when using the same 
number of time steps with the ODE45 (248 steps). However, 
when the number of steps was increased by a factor of four, i.e. 
from 248 to 992, then the time response using RK4 becomes 
more reasonable as shown in Figure 4. Concerning RK2, 



similar accuracy is obtained when the number of steps increases 
by a factor of 16 to 3968 (= 16><248 = 4x992). 
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Figure 4. Calculated residuals (Ri,R2,R3,Ra) based on Eq(51), using 
RK4 (992 equal steps) for the Lagrange. 

Also, the Crank-Nicolson method (Appendix D) based on 
two slightly different schemes leads to similar results than the 
RK4 scheme for the same number of time steps. For both 
schemes (Scheme-1 and Scheme-2), the minimum multiple of 
248 to obtain acceptable results is 4x248=992 steps (very 
similar to RK4). The results are shown in Figure 5 and Figure 
6 for scheme-1 and scheme-2, respectively. 
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Figure 5. Calculated residuals (i? 2,^3,^4) based on Eq(51), using 
Crank-Nicolson Scheme-1 for the Lagrange formulation (992 equal 
steps). 
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Figure 6. Calculated residuals (R [,R 2 ,Ri,R4) based on Eq(51), using 
Crank-Nicolson Scheme-2 for the Lagrange formulation (992 equal 
steps). 

In the sequence, in order to include the constraint of constant 
total energy, we apply the two scenarios mentioned in Section 
4.3. The results are as follows. 

Scenario-1: As an initial state we considered the numerical 
solution at the time t = T„l% (T„ = nutation period), as obtained 
using the ODE45 MATLAB implementation of the Lagrange 
formulation for very small tolerance (RelTol=le-14, 
AbsTol=le-14). It is noted that all four residuals were of the 
order 10" 14 but numerical solution crashed after T„/8 at the first 
turning point. 

Scenario-2: Now the initial conditions were taken again equal 

to <f> =0,^ =0, 6 = tt/\6,6 = 0,and y/ = ;z/l6,y> =1. 

For a better control of the computer program, the latter was 
developed in conjunction with the in house RK4. No problem 
was noticed and all residuals were found at the limit of 
computer accuracy as shown in Figure 7. 
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Figure 7. Calculated residuals (R i,i?2,^3,^4)based on Eq(51), using 
RK4 (992 equal steps) for the Lagrange formulation (Scenario-2: 
Eq(42)) in conjunction with a-posteriori correction for energy 
conservation using Eq(46). 

The time history for the three Euler angles </),0,y/ as well 
as for the level z. c of the center of mass is shown in the typical 
Figure 8; no visually significant differences were noticed for 
all models. It is worth-mentioning that the numerical solution 
predicts the extreme (turning) point at 1 74.8, which is identical 
to the aforementioned analytical value mentioned in Section 
6.1 .1 [Eq(27g) or Eq(27m)]. It is again noted thatEq(29c) is not 
applicable. 
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Figure 8. Euler angles and distance of the center of mass from the 
horizontal using ODE45 MATLAB solution with 248 steps of variable 
size. 

Moreover, the time history of the components of the support 
forces are shown in Figure 9 where one can notice that the 
vertical is initially upward (positive) until the body z axis 
reaches the horizontal and then it becomes downward (negative) 
until the lowest turning point. In the bottom right of the same 
figure, we also include the total energy and its distribution in 
terms of potential and kinetic components. 
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Figure 9. Support forces and energy using ODE45 MATLAB solution 
with 248 steps of variable size. 

As previously mentioned, the difference between the 
analytical [Eq(29c) until Eq(32b)] and the numerical solution 
was found to be enormous. In other words, the aforementioned 
analytical solutions are not applicable to a slow top as the case 
of this example. 

6.2 Example 2: The spinning wheel 

The meaning of this example is that it offers the possibility to 
modify parameters in the real world. It refers to a spinning 
flywheel, which is mounted at the end of a massless rigid rod of 
length / as shown in Figure 10. The flywheel is a cylindrical 
body of radius r and thickness h. The rigid rod is connected at 
the center of mass of the aforementioned cylindrical body. 
Under these circumstances, the momenta of inertia with respect 
to the origin O are given by: 



/ = — (3r +h ) + ml , L = —r 
12 v ' 2 



We choose the following data: 

Length of rod: / = 1 .0 m 

Gravitational acceleration: g = 9.81 m/sec 2 

Thickness of spinning wheel: h = 0.010 m 

Density: p = 7800 kg/m 3 

Radius of spinning wheel: r = 0.10m (see Fig. 10) 



(52) 



Initial conditions: (j) Q = 0,^ = , O = 7v/16,0 o = , 
and y/ =0,y/ = variable . 



gravity 



spinning wheel 




Figure 10. Spinning wheel. 

According to Eq(28), the characterization of a "fast top" 
occurs when J (a 3 » 4^ , or equivalently [by virtue of (52)]. 

Therefore: 

Condition for fast top : a - » ggf t r 2 . (53) 

For the particular data, the critical angular velocity of the 

2 i 
spinning wheel equals to C0 3 CR = 7848 s (therefore 

a >>fiV s =SS.6 s 1 ^ 846 RPM). 

In order to validate the analytical solutions, numerical 
solutions are derived for several multiples and submultiples of 



the 



critical 



angular 



velocity 



co, 



Mass of wheel: 



-phnr ■ .(= 2-45 kg) 



3, CR 

( <y 3 = Ag) 3 CR , 0. 1 < A < 10 ) using ODE45 MATLAB (for 

the default) in conjunction with Lagrange formulation (Section 
4.2). Table 1 presents the comparison of results obtained using 
analytical and numerical solution. Therefore, if the average 

angular velocity (j) is estimated by Eq(32b), the calculated 
period does not always correspond to azimuthal angle of 360 
degrees. For example, one can notice in Table 1 that at the 

critical condition ( 0) 3 =C0 2CR ) it corresponds to 428.1 degrees, 

at a double value to 373.1 degrees, whereas at a half value to 
338.2 degrees. It can be also noticed that when the spin velocity 

is smaller than C0 JCR l2 the center of mass may travel below 

the horizontal level, whereas when C0 3 > C0 3 CR the nutation is 
small and for still higher spins it tends to vanish. 
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Table 1 : Evaluation of analytical solutions 



Multiple or 
sub -multiple of 

(A) 


Angular velocity 


Maximum angle 
(^) in degrees, 
based on (32b) 


Minimum 
Level 

S*sm 


Maximum 

Level 

ferns 




(KPM) 


3.1 


8.86 


84.6 


6.9 


-0.9606 


0.9808 


0.2 


17.72 


169.2 


76.9 


-0.S425 


0.9808 


0.3 


26.58 


253.8 


229.5 


-0.6450 


0.9808 


0.4 


35.44 


338.4 


293.5 


-0.3692 


0.9808 


0.5 


44.29 


423.0 


338.2 


-0.0224 


0.9808 


0.6 


53.15 


507.6 


410.0 


0.3899 


0.9808 


0.7 


62.01 


592.2 


503.8 


0.7822 


0.9808 


0.8 


70.87 


676.8 


484.3 


0.9220 


0.9808 


0.9 


79.73 


761.4 


428.8 


0.9513 


0.9808 


1.0 


88.59 


846.0 


428.1 


0.9621 


0.9808 


1.1 


97.45 


930.6 


404.1 


0.9675 


0.9808 


1.5 


132.88 


1268.9 


384.8 


0.9755 


0.9808 


2.0 


177.18 


1691.9 


373.1 


0.9781 


0.9808 


3.0 


265.77 


2537.9 


365.8 


0.9797 


0.9808 


4.0 


354.36 


3383.8 


363.1 


0.9802 


0.9808 


5.0 


442.94 


4229.8 


361.6 


0.9804 


0.9808 


6.0 


531.53 


5075.8 


361.1 


0.9805 


0.9808 


10.0 


885.89 


8459.6 


360.5 


0.9807 


0.9808 


20.0 


1771.78 


16919.2 


360.1 


0.9808 


0.9808 



Eq(29c) up to Eq(32b). Basically, there are four physical 
quantities to calculate: 

1) The change A6 of the inclination angle during 
nutation. Equivalently, the range in which the 
z-coordinate of the centroid varies. To this, Eq(29c) is 
applied. 

2) The frequency k of the nutation (Eq(30)). 

3) The average frequency (equivalently, of angular velocity) 
of the precession (Eq(32b)). 

4) The average number of nutations in a whole period of 
precession (360 degrees of <f>). The latter equals to the 
ratio of the average frequency of precession over the 
frequency of nutation. Therefore, by virtue of Eq(30) and 
Eq(32b), the number of nutations in a whole period of 
precession are approximated by: 

Number of nutations per period of precession: 

N_ ,_. = v " (54) 



I x mgl 

For example, in case ofl=2, Eq(29) gives U 2 —U l « 0.0024 , 

which is very near to the value 0.0027 that was calculated by 
the ODE45. Also, Eq(54) predicts 31.9 nutations (30.8 when 
using the correcting factor 360/373.1 according to Table 1), 
which is very close to the measured 29 nutations as shown in 
Figure 11. It is noted than calculations concern 373.1 degrees 
of precession, whereas the overlapping of one period of 
nutation can be noticed near the initial point which is located at 

(W 0! z )« (0,-0.195, 0.981). 




Figure 11. Calculated nutations when the initial spin angular velocity 
is twice larger than the critical value of 88.59 s" 1 . 

7. DISCUSSION 

This study focuses on the "spinning top", which works 
equivalently with a flywheel at the end of a rigid rod. The 
common feature is that in both cases the distance of the centroid 
from the support is different than zero (1^0). Although in 
engineering praxis the so-called "gyroscope" is of higher 
importance, the only essential difference between a spinning 
top and a gyroscope lies on the fact that a gyroscope is usually 
supported at its centroid (thus / = ). 

It was shown that, in principle, the ODE45 MATLAB 
function can be safely applied for the numerical integration of 
the differential equations of motion, in all formulations. The 
advantage of this function is that it continuously varies the time 
step so as to locally achieve the given tolerance. 

In general, basic conclusions can be derived on the basis of 
simple closed form analytical solutions provided we carefully 
deal with the slow top and fearless with a fast top in which the 
initial spin is higher than its critical value. In the case of a 
spinning flywheel of cylindrical shape, a simple formula was 
given for the aforementioned critical value in Eq(53). It can be 
noticed than in the latter case the critical spin depends only on 
the length / and the radius r and not on the thickness h of the 
cylindrical wheel. 

8. CONCLUSIONS 

In this study we successfully applied several numerical 
schemes for the solution of the equations of motion for a 
symmetrical spinning top, which was considered as a rigid 
body with a fixed and frictionless point. It was found that using 
the standard MATLAB ODE45 with default accuracy, the 
Lagrange based formulation requires a smaller number of time 
steps than the Euler equations formulation (about 36 percent 
reduction). It was also found that the same occurs when 
increasing the tolerance variables. Both of the tested 
Crank -Nicolson schemes did not perform better than in house 
Runge-Kutta RK4 when applied in conjunction with an 
invariable size of the time step. Although the ordinary 
differential equation produced by the energy conservation is 
solvable in terms of the inclination Euler angle, it was found to 
be trapped between two "turning points". In contrast, the 
accuracy was significantly increased when the Lagrange 
solution was combined with a-posteriori correction of the time 



derivative of the inclination angle for each time step. It was 
clearly shown that although some analytical solutions can be 
applied to both slow and fast tops, the behavior of fast spinning 
tops is predicted by much simpler formulas. 
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APPENDIX A 

Derivation of motion Euler equations using matrices 

In general, the relation between the co-ordinates of a vector v in 
the local (body) Oxyz system and the global (space) Oxiy.^x 
system, is given by the simple formula: 



' local 



.(A- 



Y global ) ^\ global + A ' ^global 



and 



V global ( A ■ V local ) ( A ) V local + A ^ local ' 



By virtue of (A-5), Eq(A-7) becomes; 
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(A-6) 



(A-7) 



V local A ' y global ■ 



(A-l) 



v = A' 

V global n - 



' ^local ' V local + A "^/oca/ ' 



(A-8) 



where A is the transformation matrix from the space (global) to 
the body (local) system. In our case, the matrix A is given by 
Eq(5). 

Equivalently, it holds: 



global 



local ' 



(A-2) 



where A is the transpose of the matrix A (obviously, 

AA r =AA-'=I). 

The local Cartesian system may be considered that is 
produced by the global through a rotation G = \ e 6 around the 

unit vector \ g . Since the rotation vector 9 can be a function 

of time, i.e. 8 = 9 (t ) , then the angular velocity G) of the rigid 
body (local system) will be: 



(9 



cl() 

dt 



My 

CO, 



(A-3) 



If we redefine the abovementioned angular velocity in a 
tensor from: 





' 


-co z 






^local - 


co z 







(A-4) 




-co y 




; 





it can be easily proved that 
A = Q? local • A , or equivalently, ( A) = A T ■ CI 



local 

(A-5) 



It is worth-mentioning that the definition of the matrix £2 in 
(A-4) has been chosen so as it fulfills the identity: 

fiv = wxv (matrix product equals to the cross product). 

Therefore, the components of the time derivative of the 
vector v are related by: 



Let us now designate by L the angular momentum of the 
rigid body and M the external moment (or torque). The second 
Newton's law for the rotation is written as follows: 



iyi _ ^global _ j 

Lyl global ^ - '-'global 



(A-9) 



Since both M and L are vectors, they obey the general 
transformation which is given by Eq(A-l) and (A-6). Therefore 
it holds: 



M = A -M 

global A iTA local' 



^global A- 



(A-10) 
(A-ll) 

Substituting (A-10) and (A-ll) into (A-9), and considering 
the identity (A-7), one obtains: 

A r -M local = (A 7 )* L local + A T t local (A-12) 
Substituting (A-5) into (A-12) one obtains: 

A r • M ,oca, = A 7 • ^loca, ' ^local + ^ ' ^ local 

(A- 13) 

Left-multiplying (A-13) by A , considering that AA 7 " = I , 
and rearranging the left and right members, we receive: 

Kcal + ^local ■ ^local = M ,oca, ( A "14) 

Considering the notation appearing in the main text, i.e.: 



(A- 15) 





w 






L = 


I 2 co 2 


, and L = 


I 2 C0 2 




I 3 co 3 




I 3 C0 3 



^local - 



CO, 



K -co 2 



-co } 


CO, 



CO 



(A- 16) 



and 



M 



local 



M 2 



(A- 17) 



Equation (A-14) finally becomes: 



/,<», 




I 2 cb 2 


+ 


I 3 cb 3 


V 



CO, 



-co 3 co 2 
-co, 



-co 2 co l 







\ 


l x co x 




X" 




I 2 co 2 




M 2 


/ 


J.co, 







y„ +1 =y„+k 2 +o(^ 3 ) 
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(C-3) 



(A-18) 



k 2 =hf 



h k, A 



Equation (A-18) is identical with Eq(3a,b,c), and this 
completes the proof. 

APPENDIX B 

The angular velocity in the fixed co-ordinate 

system 

Applying (A-2) for the vector \ local = \co { G) 2 G) 3 ] , of 
which the components in the fixed system is 

V g i oba , = CO yl CO zl ] , we finally obtain: 



co xl = 6 cos cf> + y/ sin cf> sin , 

co vl =6 sincf>-y/ coscf>sin6 , (B-l) 

co zl =cj) + y/ cosO . 

It is remarkable that the same result will be obtained if we 
transform the tensor given in (A-16) to the Q g i oba i, 

using the formula: 



global 



global 



(B-2) 



Also, the four-point Runge-Kutta (RK4) solution is produced 
using the recursion: 

" 6 3 3 6 K ' 
K-hf{t n ,y n ), 

( h k "\ 
k,=M (+-y ii+ ^ , (C-4) 
v 2 2 J 

k 4 = Af (f. + A.y. + k^) 
Both (C-3) and (C-4) restrict to the same (constant) time step h. 
APPENDIX D 
Crank-Nicolson method as applied 

For the problem described by (C-l) and (C-2), the trapezoidal 
rule suggests using: 



APPENDED C 
Runge-Kutta algorithms as applied 

For a first order differential system 

y=f(*,y), (c-i) 

with initial conditions 

y(0) = y , (C-2) 



the two-point Runge-Kutta (RK2) solution is produced using 
the recursion: 



y n+P =(i-/3)y„+/?y„ + p 
y n+P =(y n+ i-y„)/ A ' 



(D-l) 



and 



fjt,y) = (l-p)f n +pf n+i (D-2) 
Therefore, the system (C-l) is solved using the recursion: 

y„ + i=y.+A[(i-£)f.+/». + J (°- 3 ) 

According to the choice of the parameter J3 we distinguish 
four alternative schemes as follows: 



P = 



0, forward difference / Euler 
1/2 , Crank-Nicolson (mid-point) 
2/3, Galerkin 

1, backward difference 



(D-4) 



In this paper we have tested only the case of /? = 1/2 , which is 
the so-called Crank-Nicolson scheme. Moreover, since the 



right-hand-side f (^y) is a function of the unknown value y, 
we have to start (D-3) for the initial given value f„ at t=0, and 
then perform some iterations to update f , which is a 

function of y jl+1 . We distinguish two alternative schemes: 

1) Scheme 1: We use the formula: 

and perform a number of iteration to determine the 
right-hand-side term {f n + f n+] ^ J 2 . 

2) Scheme 2: As an alternative, we use the modified 

formula: 



(D-6) 



and perform a number of iteration to determine the 



( 



right-hand-side term f 
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